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1.  INTRODUCTION 


3.  FOURTH-ORDER  DISCRETIZATION 


Among  the  many  sources  of  error  in  modeling  tropical 
cyclones  is  discretization  error:  how  well  are  the  governing 
equations  approximated  numerically?  While  discretization 
error  is  typically  small,  numerical  methods  with  smaller  er¬ 
rors  should  be  preferred  over  others,  all  else  being  equal. 
This  paper  addresses  the  tradeoff  between  higher-accuracy 
discretization  and  increased  computational  work  in  an  adap¬ 
tive  multigrid  tropical  cyclone  track  model. 

2.  MODEL  DESCRIPTION 

The  model  described  here  is  a  higher-order  version  of 
the  MUDBAR  model  described  in  Fulton  (1997,  2000). 
Formulated  on  a  section  of  the  sphere  using  a  Mercator 
projection  (true  at  latitude  <p  =  (pc),  the  model  consists  of 
the  modified  barotropic  vorticity  equation 

+m2J(V',C) =  i/m2V2C,  (1) 

with  relative  vorticity  £  and  streamfunction  ip  related  by 


Here  V2  =  d2 /dx2  H-  d2/dy2,  J(V>,0  is  the  JacobiarT 
of  (ipX)  with  respect  to  (xyy),  /?  =  20  cos  (p/a  (with  a 
and  O  the  radius  and  rotation  rate  of  the  earth),  m  = 
cos  <pd  cos  <p  is  the  map  factor,  v  is  the  diffusion  coeffi¬ 
cient,  and  7  is  the  inverse  of  the  effective  Rossby  radius. 
The  model  domain  is  a  rectangle  in  x  and  y  centered  at 
(x,y)  —  (0,0),  where  (A,  <p)  =  (Ac,<£c).  At  the  bound¬ 
aries  we  specify  the  streamfunction  ip  (and  thus  the  normal 
component  of  the  velocity);  where  there  is  inflow,  we  also 
specify  the  vorticity  £. 

The  equations  are  discretized  in  space  by  finite  dif¬ 
ferences  on  uniform  rectangular  grids,  approximating  the 
advection  terms  by  the  Arakawa  Jacobian.  The  model 
achieves  higher  resolution  near  the  vortex  by  superimposing 
nested  overlapping  grids  with  different  mesh  sizes.  Unlike 
conventional  nested-grid  methods,  using  multigrid  process¬ 
ing  (Brandt,  1977)  in  solving  for  the  streamfunction  allows 
optimal  solution  speed  and  accurate  estimates  of  trunca¬ 
tion  error;  the  latter  can  be  used  in  an  adaptive  mesh  re¬ 
finement  scheme  to  provide  just  the  resolution  needed  at 
each  point.  More  detail  for  the  second-order  model  is  given 
in  Fulton  (2000);  here  we  investigate  the  gains  possible  by 
using  fourth-order  differencing. 
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The  model  is  discretized  in  time  using  the  classical 
fourth-order  Runge-Kutta  scheme  (also  used  in  the  second- 
order  model).  To  achieve  fourth-order  discretization  in 
space  we  use  the  methods  described  below. 


3.1  Streamfunction  equation  solution 


For  the  Helmholtz  problem  (2)  we  use  a  compact 
fourth-order  discretization  known  as  the  Mehrstellen  Ver- 
fahren  (MV)  discretization  (e.g.,  Schaffer,  1984).  For  the 
Poisson  problem  (i.e.,  setting  7  =  0  and  m  -  1)  this  takes 
the  form 
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Here  iph  and  £h  denote  grid  functions  consisting  of  the  ap¬ 
proximate  values  of  streamfunction  and  vorticity  on  a  uni¬ 
form  rectangular  grid  with  mesh  spacing  h  in  x  and  y,  and 
the  square  brackets  denote  difference  stencils  on  that  grid. 
Since  this  compact  discretization  involves  only  the  nearest- 
neighbor  points,  it  can  be  applied  at  all  interior  points  on 
the  grid  (with  Dirichlet  boundary  conditions  specifying  the 
boundary  values). 

To  solve  the  resulting  discrete  problem  efficiently  we 
use  a  multigrid  method.  As  in  the  second-order  model,  this 
uses  the  Full  Multigrid  (FMG)  approach,  solving  on  each 
grid  level  from  coarsest  to  finest  in  turn,  interpolating  the 
solution  from  the  previous  level  to  provide  the  initial  approx¬ 
imation  on  the  new  level.  The  fourth-order  discretization 
(3)  is  used  only  on  the  “currently  finest”  level.  Compared 
to  the  second-order  method,  three  changes  are  needed: 

(1)  Gauss-Seidel  relaxation  is  used  with  lexicographic  or¬ 
dering  (rather  than  red-black); 

(2)  Two  V(2,l)  cycles  are  needed  to  solve  on  each  level 
(rather  than  one  V(l,l)  cycle); 

(3)  The  initial  interpolation  to  a  new  level  requires  sixth- 
order  accuracy  (rather  than  fourth-order  bicubic  inter¬ 
polation). 

The  change  in  relaxation  ordering  is  a  result  of  using  a  nine- 
point  stencil;  the  other  two  changes  are  required  to  ensure 
that  the  discrete  problem  is  solved  adequately  on  any  given 
grid  and  that  its  solution  is  accurately  approximated  on  the 
next  grid.  With  these  changes  we  find  that  the  resulting 
method  [in  multigrid  terminology,  a  2-FMV(2,l)  method] 
solves  to  the  level  of  truncation  error  (i.e.,  residual  norm  less 
than  the  truncation  error  norm)  for  problems  with  known 
(smooth)  solutions,  independent  of  the  mesh  size  h.  Nu¬ 
merical  tests  verify  that  the  error  in  the  computed  solution 
of  the  Helmholtz  problem  is  indeed  0(h 4). 
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3.2  Vorticity  equation  discretization 

The  standard  second-order  Jacobian  of  Arakawa  (1966) 
uses  only  nearest-neighbor  points  and  conserves  discrete 
analogues  of  vorticity,  enstrophy,  and  kinetic  energy.  Un¬ 
fortunately,  it  can  be  shown  that  no  fourth-order  Jacobian 
has  the  same  properties:  to  achieve  fourth-order  accuracy 
while  maintaining  the  conservation  properties  one  must  in¬ 
clude  points  other  than  the  nearest  neighbors.  Thus,  we  use 
the  fourth-order  Jacobian  from  the  Appendix  of  Arakawa 
(1966)  at  all  interior  points  on  the  grid  except  those  along 
the  first  interior  lines  (adjacent  to  each  boundary),  where 
we  use  the  usual  second-order  Jacobian.  Since  the  set  of 
points  where  lower  accuracy  is  used  is  small,  this  has  neg¬ 
ligible  effect  on  the  overall  accuracy.  Numerical  tests  with 
specified  ip  and  £  verify  that  the  truncation  error  for  the 
Jacobian  is  indeed  0(h 4). 


Does  the  asymptotic  superiority  of  fourth-order  differ¬ 
encing  pay  in  practice?  To  answer  this,  we  ran  both  models 
with  variety  of  grids  (uniform  grids  and  local  refinements 
using  patches  of  various  specified  sizes).  Figure  2  shows  the 
resulting  mean  track  errors  as  a  function  of  the  computer 
time  (for  a  72  hour  model  run).  For  the  uniform-grid  runs 
(points  joined  by  lines  in  Fig.  2)  the  fourth-order  model 
is  approximately  5-10  times  faster  than  the  second-order 
model  (for  the  same  level  of  error);  for  the  same  amount  of 
computer  time,  it  is  at  least  5  times  more  accurate.  With 
local  refinement  (remaining  points  in  Fig.  2),  fourth-order 
differencing  shows  some  advantage  for  moderate  resolution 
(e.g.,  errors  of  10-30  km);  however,  the  advantage  appar¬ 
ently  disappears  at  higher  resolutions.  This  latter  result 
may  reflect  the  need  for  higher-order  grid  transfers  in  the 
fourth-order  model  with  local  refinement. 


4.  NUMERICAL  RESULTS 

To  quantify  the  model  performance  we  use  the  track 
error,  defined  as  the  difference  in  vortex  position  compared 
to  a  high-resolution  reference  run  ( h  =  4  km).  Figure  1 
shows  the  track  error  as  a  function  of  h  for  several  uniform- 
grid  runs  of  the  second-order  and  fourth-order  models  at  24, 
48,  and  72  hours  and  the  mean  errors  over  each  72  hour 
model  run.  The  slopes  of  the  curves  verify  that  the  overall 
model  converges  at  approximately  the  proper  rate. 


Figure  1.  Overall  convergence  of  the  model: 

(a)  second-order  version,  (b)  fourth-order  version. 


Figure  2.  Track  error  vs.  CPU  time. 
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